Table of Contents¶
import IPython.display as Display
from IPython.display import Latex
%matplotlib inline
%run ../starter.py
from PIL import Image
from scipy import signal as sig
import matplotlib.lines as mlines
from matplotlib.colors import LogNorm
Escletxes verticals¶
Estudi de les interferències fetes per una - o mútliples - escletxes verticals.
Extracció de les dades¶
Carreguem les dades¶
Hem fet fotografies als papers milimetrats on hem pres les dades. A partir de les fotografies captarem els mà xims, i podrem fer les mesures pertinents. Començem carregant les coordenades en la imatge dels papers milimetrats de les dades preses:
data = pd.read_csv("IN/CSV/positions.dat", sep="\t")
data["edges"] = [((X0,Y0),(X1,Y1))
for X0,Y0,X1,Y1 in
zip(data['X0'],data['Y0'],
data['X1'],data['Y1'],)]
data
no_secondary_images = data[data["N"]==1]
secondary_images = data[data["N"]>1]
Imatges preses que estudiarem:
from io import BytesIO
for im, edges in zip(data['IMG'],data['edges']):
im = load_subImage(im, edges,"IN/IMG/im0{}.jpg")
im = Image.fromarray(im)
bio = BytesIO()
im.save(bio, format='png')
Display.display(
Display.Image(bio.getvalue(), format='png'))
Ara que tenim les coordenades, podem construir una funció que ens permeti llegir les imatges, i extreure'n les dades:
def load_subImage(imNum, coords,imName = "IN/IMG/im0{}.jpg"):
img = Image.open(imName.format(imNum))
img = np.asarray(img)
return img[coords[0][1]:coords[1][1],
coords[0][0]:coords[1][0],:]
Un exemple de les imatges fetes podria ser la següent:
.
Ara ens interessaria calcular de forma automà tica la posició dels mà xims en les imatges.
Càlcul dels màxims en les imatges¶
Com es pot observar, les imatges tenen molt de soroll, per tant cal filtrar les imatges (qualsevol programa d'edició d'imatges de propòsit general aniria bé). Un cop correctament filtrada, la imatge es transforma en:
.
Per a trobar les barres de la figura, utilitzarem el següent algoritme:
Projectar la imatge sobre l'eix de les X: Sumar totes les columnes de pÃxels en l'eix Y. Aixà obtindrem una tira 1-dimensional de nombres. El valor d'aquesta tira ens ha de permetre detectar els punts més foscos de la imatge.
Filtrar la sèrie obtinguda de dades, ja que tindrà molt soroll.
Detectar els mà xims obtinguts després del filtre
La següent funció implementa aquest algoritme:
def find_maximum(img, N=5, M=35, pc=75):
#Projecció
lkh = (255-img[:,:,2]).sum(axis=0)
#filtre
smooth1 = np.convolve(lkh, np.ones((N,))/N, mode='same')
smooth2 = sig.medfilt(lkh, M)
smooth =(smooth1-smooth2)[25:-25]
#detecció de mà xims
sel1 = (np.diff(smooth,1)[:-1]>0)&(np.diff(smooth,1)[1:]<=0)
sel2 = np.percentile(smooth, pc) < smooth[1:-1]
mx = sel1&sel2
return mx, smooth
Provem per tant, d'aplicar aquest algoritme a les dades donades. Posem dos exemples dels resultats, un per al cas en que hem marcat els zeros de la ona moduladora, i un per al cas en que no hi havia ona moduladora:
im = load_subImage(secondary_images['IMG'][8],
secondary_images['edges'][8])
mx, smth = find_maximum(im)
rg = np.array(range(len(smth)))
plt.plot(smth, lw=.5, c='black')
plt.plot(rg[1:-1][mx], smth[1:-1][mx],'.', color='#880000', ms=10)
plt.title(u"Projecció de les intensitats de la imatge\n"+
u"Cas amb ona moduladora\n", fontsize = 15)
plt.xlabel(u"Posició (Px)")
plt.ylabel(u"Projecció de les intensitats")
plt.ylim([-100,8000])
plt.xlim(0,1500)
plt.subplots_adjust(top=0.85)
plt.savefig("OUT/FIG/point_detector_mult.png", dpi = 300)
im = load_subImage(no_secondary_images['IMG'][2],
no_secondary_images['edges'][2])
mx, smth = find_maximum(im)
rg = np.array(range(len(smth)))
plt.plot(smth, lw=.5, c='black')
plt.plot(rg[1:-1][mx], smth[1:-1][mx],'.', color='#880000', ms=10)
plt.title(u"Projecció de les intensitats de la imatge\n"+
u"Cas sense ona moduladora\n", fontsize = 15)
plt.xlabel(u"Posició (Px)")
plt.ylabel(u"Projecció de les intensitats")
plt.ylim([-100,4000])
plt.xlim(0,1500)
plt.subplots_adjust(top=0.85)
plt.savefig("OUT/FIG/point_detector_sing.png", dpi =300)
Ara que tenim la detecció dels mà xims feta, queda una petita feina a fer per al cas en què hi hagi ona moduladora: Discriminar els punts referents a la ona moduladora dels referits a la ona modulada. Cal dissenyar, per tant, un algorisme per a fer-ho.
Separació dels màxims en les imartges¶
Per a fer-ho, tal i com ens fa pensar la intuició, senzillament posarem un treshold, i agafarem les que siguin superiors al treshold.
La següent funció implementa aquest algoritme:
def separate_maximums(smth,mx, q=0.5):
rmx = smth[1:-1][mx]
m1 = np.sort(rmx)[-3]
m2 = np.percentile(rmx, 75)
treshold = m1*q+m2*(1-q)
#treshold = m2
return np.logical_and(mx,smth[1:-1]>treshold), np.logical_and(mx,smth[1:-1]<=treshold)
Mostrem ara el resultat de l'algorisme implementat per al mateix cas que estavem treballant:
im = load_subImage(secondary_images['IMG'][8],
secondary_images['edges'][8])
mx, smth = find_maximum(im)
mx1, mx2 = separate_maximums(smth,mx)
rg = np.array(range(len(smth)))
plt.plot(smth, lw=.5, c='black')
plt.plot(rg[1:-1][mx1], smth[1:-1][mx1],'.', color='#880000', ms=10)
plt.plot(rg[1:-1][mx2], smth[1:-1][mx2],'.', color='#008800', ms=10)
plt.title(u"Projecció de les intensitats de la imatge\n"+
u"Cas amb ona moduladora\n", fontsize = 15)
plt.xlabel(u"Posició (Px)")
plt.ylabel(u"Projecció de les intensitats")
plt.ylim([-100,8000])
plt.xlim(0,1500)
plt.subplots_adjust(top=0.85)
plt.savefig("OUT/FIG/point_detector_separator.png", dpi = 300)
Càcul de les distàncies entre pics (en píxels)¶
Ara ens hem d'adonar d'un punt molt important a tenir en compte: les dades que hem trobat no son, ni de bon tros, perfectes. Per tant, no podem calcular la distà ncia entre mà xims fent diferència entre mà xims consecutius, ja que el resultat fallarà bastant estrepitosament allà on tinguem mà xims que sobrin/faltin. Per tant, haurem de reomplir els forats que faltin. Per a fer-ho, detectarem on falten punts, i assignarem a cada punt un index $k$, que ens digui quin punt seria si els tinguessim tots:
def find_real_point_position(positions):
orders = range(len(positions))
m = np.mean(np.diff(positions)/np.diff(orders))
orders = np.cumsum(np.round(np.diff(positions)/m))
orders = np.append((0),orders)
return orders
for im,edge in zip(secondary_images['IMG'],
secondary_images['edges'],):
im = load_subImage(im,edge)
mx, smth = find_maximum(im)
mx1, mx2 = separate_maximums(smth,mx)
rg = np.array(range(len(smth)))
points_zero = rg[1:-1][mx2]
plt.plot(points_zero-points_zero[0],'o')
plt.title(u"Posició dels zeros de la ona moduladora\n"+
u" en funció del nombre del"+
u" mà xim (original)\n", fontsize=15)
plt.xlabel(u"#mà xim")
plt.ylabel(u"Posició (px)")
plt.subplots_adjust(top=0.85)
plt.savefig("OUT/FIG/zeros_before_filter.png", dpi = 300)
for im,edge in zip(secondary_images['IMG'],
secondary_images['edges'],):
im = load_subImage(im,edge)
mx, smth = find_maximum(im)
mx1, mx2 = separate_maximums(smth,mx)
rg = np.array(range(len(smth)))
points_zero = rg[1:-1][mx1]
plt.plot(points_zero-points_zero[0],'o')
plt.title(u"Posició dels zeros de la ona moduladora\n"+
u" en funció del nombre del"+
u" mà xim (original)\n", fontsize=15)
plt.xlabel(u"#mà xim")
plt.ylabel(u"Posició (px)")
for im,edge in zip(secondary_images['IMG'],
secondary_images['edges'],):
im = load_subImage(im,edge)
mx, smth = find_maximum(im)
mx1, mx2 = separate_maximums(smth,mx)
rg = np.array(range(len(smth)))
points_zero = rg[1:-1][mx2]
positions = find_real_point_position(points_zero)
plt.plot(positions,
points_zero-points_zero[0],'o')
plt.title(u"Posició dels zeros de la ona moduladora\n"+
u" en funció del nombre del"+
u" mà xim (corregit)\n", fontsize=15)
plt.xlabel(u"#mà xim")
plt.ylabel(u"Posició (px)")
plt.subplots_adjust(top=0.85)
plt.savefig("OUT/FIG/zeros_after_filter.png", dpi = 300)
for im,edge in zip(secondary_images['IMG'],
secondary_images['edges'],):
im = load_subImage(im,edge)
mx, smth = find_maximum(im)
mx1, mx2 = separate_maximums(smth,mx)
rg = np.array(range(len(smth)))
points_zero = rg[1:-1][mx1]
positions = find_real_point_position(points_zero)
plt.plot(positions,
points_zero-points_zero[0],'o')
plt.title(u"Posició dels zeros de la ona modulada\n"+
u" en funció del nombre del"+
u" mà xim (corregit)\n", fontsize=15)
plt.xlabel(u"#mà xim")
plt.ylabel(u"Posició (px)")
for im,edge in zip(no_secondary_images['IMG'],
no_secondary_images['edges'],):
im = load_subImage(im,edge)
mx, smth = find_maximum(im)
rg = np.array(range(len(smth)))
points_zero = rg[1:-1][mx]
plt.plot(points_zero-points_zero[0],'o')
plt.title(u"Posició dels zeros"+
u" en funció del nombre del"+
u" mà xim,\n sense ona modulada (original)\n", fontsize=15)
plt.xlabel(u"#mà xim")
plt.ylabel(u"Posició (px)")
for im,edge in zip(no_secondary_images['IMG'],
no_secondary_images['edges'],):
im = load_subImage(im,edge)
mx, smth = find_maximum(im)
rg = np.array(range(len(smth)))
points_zero = rg[1:-1][mx]
positions = find_real_point_position(points_zero)
plt.plot(positions,
points_zero-points_zero[0],'o')
plt.title(u"Posició dels zeros"+
u" en funció del nombre del"+
u" mà xim,\n sense ona modulada (corregit)\n", fontsize=15)
plt.xlabel(u"#mà xim")
plt.ylabel(u"Posició (px)")
Ara podem utilitzar aquesta correcció que hem trobat per a estimar la distà ncia entre mà xims adequadament. Només cal recordar que si be en mesures fÃsiques s'acostuma a utilitzar la mitjana i la desviació està ndard com a estimadors del valor i la seva possible fluctuació, en mesures de processament d'imatge no és aixÃ. Cal utilitzar esimadors més robustos, que tinguin en compte que un cert percentatge de les dades pot ser senzillament erroni, i no aportar cap informació útil. Els estimadors a utilitzar en aquests casos són la mediana i la median absolute deviation. Sabent això, procedim a calcular els estimadors dels parà metres:
modulated = unp.uarray(np.zeros(11),0)
modulator = unp.uarray(np.zeros(11),0)
for k in range(11):
im = load_subImage(data['IMG'][k], data['edges'][k])
mx, smth = find_maximum(im)
mx1, mx2 = separate_maximums(smth,mx)
rg = np.array(range(len(smth)))
if data['N'][k]>1:
mx=mx1
points_zero = rg[1:-1][mx2]
positions = find_real_point_position(points_zero)
deltaN = np.diff(positions)
deltaP = np.diff(rg[1:-1][mx2])
delta = deltaP[deltaN>0]/deltaN[deltaN>0]
med = np.median(delta)
mad = np.median(np.abs(delta-med))
modulated[k]=uc.ufloat(med,mad/np.sqrt(len(delta)))
points_zero = rg[1:-1][mx]
positions = find_real_point_position(points_zero)
deltaN = np.diff(positions)
deltaP = np.diff(rg[1:-1][mx])
delta = deltaP[deltaN>0]/deltaN[deltaN>0]
med = np.median(delta)
mad = np.median(np.abs(delta-med))
modulator[k]=uc.ufloat(med,mad/np.sqrt(len(delta)))
Els resultats que hem trobat els mostrarem més endavant, en una taula completa que ens mostri tots els cà lculs fets.
Transformació de les dades a mm¶
Ara, hem de transformar els resultats que hem trobat en pÃxels a milÃmetres. Per a fer això hem de calcular la distà ncia entre dos punts a distà ncia coneguda en la imatge. Hem escollit dos extrems de la zona milimetrada, que té $280$ mm de longitud. Amb els pÃxels trobats, les ratios associades són:
ratios = {1:(1547 - 47)/uc.ufloat(280,2),
2:(1550 - 45)/uc.ufloat(280,2)}
Ara apliquem les ratios que hem trobat a cada mesura:
ratio_list = [ratios[k] for k in data['IMG']]
modulated/=ratio_list
modulator/=ratio_list
data['modulated'] = modulated
data['modulator'] = modulator
Potser, per altra banda, ja seria hora de començar a mostrar les dades que hem obtingut! [recordem que les dades modulades són en mm]
data
Anàlisi de les dades obtingudes¶
Ara que ja tenim les dades, podem començar a comprovar les formules donades al guió. Podem calcular els angles associats als desplaçaments estudiats, tant des d'un punt de vista teòric, a partir de les mides, com a partir dels desplaçaments, per a observar si són coherents:
distance = uc.ufloat(2740,20) #distà ncia a la paret
theta_modulated = modulated/distance
theta_modulator = modulator/distance
Començarem per la ona moduladora (la generada per les escletxes i no per la interacció entre escletxes)
Ara si recordem les formules del guió, tenim que la intensitat en funció de la posició és:
\begin{equation} I(\theta) = I_0 \left ( \frac {\sin \beta}{\beta}\right ) ^2 \end{equation}on $\beta = k \frac b 2 sin \theta$
Per tant, la intensitat s'anul·la si:
\begin{equation} k \frac b 2 sin \theta = \pi n \end{equation}En aproximaciò paraxial, ho podem escriure com:
\begin{equation} \theta = \frac {2 \pi}{b k} n \end{equation}i per tant $\Delta \theta = \frac {2 \pi}{b k} = \frac{\lambda}{b}$
Si recordem que tenÃem un laser He-Ne:
lamb_HeNe = 633E-6 #en mm
plt.errorScatter( lamb_HeNe/data['Size'], theta_modulator)
plt.plot([0,0.04],[0,0.04])
plt.xlim([0,0.04])
plt.subplots_adjust(top=0.85)
plt.subplots_adjust(top=0.85)
plt.savefig("OUT/FIG/calculated_vs_nominal.png", dpi = 300)
#############^^^^^ACABAR^^^^^^##############
Anàlisi d'imatges¶
Mentre erem al laboratori và rem prendre diverses imatges, de les quals hem pogut extreure informació quantitativa, a partir de l'estudi de les imatges
Anàlisi de la forma dels pics d'intensitat¶
Un dels anà lisis més immediats és el de comprovar la forma de la intensitat dels mà xims en una difracció amb una escletxa. Per a fer-ho, hem seguit el mateix procediment que abans. Hem fet una fotografia al patró del laboratori, i ens hem restringit a la franja horizontal on es veu el patró:
img = Image.open('IN/IMG/single_slit_intensity.jpg')
img
Una vegada hem carregat aquesta imatge, seguim el procediment anterior de sumar la intensitat de totes les rectes verticals, "projectant" les intensitats sobre l'eix horizontal. Aquest resultat ara, semblaria que l'hem de representar contra la funció sinc, ja que és el que diu la formula trobada a classe. Fer-ho aixà seria desastrós.
Cal tenir en compte, per a que funcioni, un fenòmen important en les cà meres fotogrà fiques digitals: La intensitat guardada en RGBen memòria no és proporcional al nombre de fotons, si no que pateix un preprocessament important. Aquest preprocessament entre d'altres coses, acostuma a incloure una transformació logarÃtmica de les dades amb alta intensitat. Per tant, el que cal representar no és les dades contra el sinc, si no contra el seu logaritme.
Fem, per tant, el que hem comentat:
npa = np.asarray(img)[:,:,0]
k = np.sum(npa, axis=0)
x=np.linspace(-15,15,1500)
px=30./1600*np.array(range(len(k)))-14.5
plt.title(u'Intensitat en funció de la posició\n', fontsize = 15)
plt.ylabel(u'Log Intensitat\n(unitats arbitrà ries)')
plt.xlabel(u'Posició en el paper (cm)')
nk=np.convolve(k, np.ones((3,))/3, mode='same')
nk = nk/max(nk)*1.01
plt.plot(x, np.log(np.sinc(x/2.01)**2),'--', c='#880000', lw=0.5 )
plt.plot(px-.2,(nk-1.21)*10+0.03*px, c='k', label ='Valor mesurat', lw=1)
plt.ylim(-15,2)
plt.xlim(-14,14)
plt.subplots_adjust(top=0.85)
plt.savefig("OUT/FIG/intensity_from_phone.png", dpi = 300)
Anàlisi de la xarxa hexagonal¶
Más difÃcil todavÃa.... Agafarem una imatge d'un patró de difracció fet per una xarxa de difracció hexagonal, i intentarem recuperar-ne la xarxa fent la transformada de fourier inversa de la imatge. A l'esquerra mostrem la imatge, i a la dreta, la seva transformada de fourier. Seguint la filosofia d'abans, hem hagut de calcular la exponencial de la imatge, per a contrarrestar el preprocessament.
img = Image.open('IN/IMG/hex.jpg')
img
plt.figure(figsize=(8,8))
im = np.asarray(img)[:,:,0]
im=np.exp(im/25)
s=(np.log(np.abs(np.fft.fft2(im))))
plt.imshow(s[::-1], norm=LogNorm(vmin=11.5, vmax=13),extent = [0,1,0,1],cmap="Spectral")
En la imatge de la transformada de fourier es pot intuir en algunes zones una xarxa, i per les direccions de les lÃnies, podria semblar que hexagonal i tot. Ara bé, cal una mirada molt més propera al patró de difracció per a poder veure realment la existència, en alguns trossos, d'aquesta xarxa hexagonal:
def hex(r,theta):
t = 2*3.141592*np.array(range(7))/6+theta
return r*np.cos(t), r*np.sin(t)
fig1 = plt.figure(figsize=(14,16))
ax1 = fig1.add_subplot(111)
cm="Spectral"
s=(np.log(np.abs(np.fft.fft2(im))))
ax1.imshow(s[::-1], norm=LogNorm(vmin=11.5, vmax=13),extent = [0,1,0,1],cmap=cm)
plt.title(u'Transformada de Fourier del patró de difracció\n', fontsize = 18)
# add a line
ax1.plot([.4 , .4, .25,.25,.4,.6],
[0.45, .6, .6,.45,.45,.415], c='#000022', lw=2)
plt.xlim(0,1)
plt.ylim(0,1)
s=(np.log(np.abs(np.fft.fft2(im))))[450:600,250:400]
ax2 = plt.axes([.575, .185, .3, .3], axisbg='y')
ax2.imshow(s[::-1], norm=LogNorm(vmin=11.5, vmax=13),extent = [0.25,0.4,0.45,0.6,],cmap=cm)
ax2.set_xlim([0.25,0.4])
ax2.set_ylim([0.45,0.6,])
x,y = hex(0.01,0)
ax2.plot(x+0.281,y+0.518,lw=4, alpha = 1, c='k')
#ax2.plot([.3 , .4, .25,.25,.4,.6],
# [0.5, .6, .6,.45,.45,.415], c='#000022', lw=2)
for axis in ['top','bottom','left','right']:
ax2.spines[axis].set_linewidth(2)
plt.setp(ax2, xticks=[], yticks=[])
ax3 = plt.axes([.685, .642, .2, .2], axisbg='y')
ax3.set_xlabel('Imatge Original', color='w', labelpad=-20, fontsize = 15)
ax3.imshow(img)
plt.setp(ax3, xticks=[], yticks=[])
plt.subplots_adjust(top=0.85)
plt.savefig("OUT/FIG/big_fft_image.png", dpi = 300)
f = lambda x:'${:L}$'.format(x) if x.s!=0 else "--"
F = lambda x: map(f, x)
h = lambda x: str(x) if x>0 else "--"
H = lambda x: map(h, x)
theta_modulator[theta_modulator==0] = uc.ufloat(1E20,0)
theta_modulated[theta_modulated==0] = uc.ufloat(1E20,0)
r1 = pd.DataFrame({"d" : H(data["Size"]),
"sep" : H(data["Sep"]),
"modulatorl" : F(modulator),
"modulatortheta" : F(theta_modulator),
"modulatedl" : F(modulated),
"modulatedtheta" : F(theta_modulated),
"calculatedsep" : F(lamb_HeNe/theta_modulated),
"calculatedd" : F(lamb_HeNe/theta_modulator),
"N" : H(data["N"])
})
r1.to_csv("OUT/TBL/output.csv")
r1